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Abstract. We describe two separate wavelength discretization schemes that can be used in the numerical solution 
of the comoving frame radiative transfer equation. We present an improved second order discretization scheme 
and show that it leads to significantly less numerical diffusion than previous scheme. We also show that due to 
the nature of the second order term in some extreme cases it can become numerically unstable. We stabilize the 
scheme by introducing a mixed discretization scheme and present the results from several test calculations. 



1. Introduction 



The numerical solution of the radiative transfer equa- 
tion plays a large role in our interpretation of spec- 
troscopic data of astrophysical sources. New methods 
and faster computers have led to a resurgence of in- 
terest in solving the transfer equ ation (see for example 
iHubenv. Mihalas. fc Werneill2003|h 



The numerical radiative tran sfer in the co-m oving 
frame (CMF) method discussed in lHauschildtl l|l992f) uses 
a discretization of the dXI/dX terms in the RTE to obtain 
a formal solution. We show that a second order discretiza- 
tion scheme for the wavelength derivative leads to better 
numerical accuracy for a number of applications, such as 
stellar winds. In addition, the new discretization allows us 
to include the effects of the wavelength derivative in the 
construction of the approximate A operator used in the op- 
erator splitting method. This improves the computational 
performance of the algorithm. It is possible to mix the two 
discretization scheme to tailor the performance of the al- 
gorithm to the problem being considered. In the following 
we describe the new discretization and the construction 
of the A* matrix for arbitrary bandwidths and discuss 
some test results. We have impleme nted this scheme into 
the PHOENIX code (see for example IHauschildt fc Baronl 
Il999j) . 



2. Method 

In the following dis cussion we use the same notation as 
in lHauschildtl l)1992j) and reproduce the key equations for 
convenience. We use the spherically symmetric form of 
the special relativistic, time independent (d/dt = 0) RTE 
(hereafter, SSRTE), the restriction to plane parallel ge- 
ometry is straightforward. The calculation o f the charac- 
teristics is identical that of lHauschildtl l|l992j) and we thus 
assume that the characteristics are known. First, we will 
describe the process for the formal solution, then we will 
describe how we construct the approximate A operator, 
A*. 



2.1. Radiative Transfer Equation 

In t he CMF the SSRTE in the wavelength ( A) scale is gi ven 
by llMihalas fc Weibel-Mihala3ll984lHauschildlJ ll992^ 



dl 

dr 



dl 



d\I 



a r— + + a >^TT + 4a A^ = V-Xl 



with 

a,, = j(fi + /3), 
«p = 7(1 - A* 2 ) 
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r or 



a\ = 7 
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Along the characteristics, Eq. ^ has the form jMihalad 
1980) 



absorptive continuum, no 



dh 
ds 



9X1 , A ST 

a l~g^ = 7 ll - \Xl + 4a/)I; 



(5) 



where ds is a line element along a (curved) characteris- 
tic, Ii(s) is the specific intensity along the characteristic 
at point s > (s = denotes the beginning of the char- 
acteristic) and wavelength point A;. The coefficient a/ is 
defined by 



a; = 7 



/3(1-M 2 ) 



7V (A* + 0) 



dp 
dr 



where = v/c,j = y/l — P 2 and r is the radius. i]i and xi 
are the emission and extinction coefficients at wavelength 
A/, respectively. 

2.1.1. First dl/dX discretization 

As in iMihalas fc Weibel-Mihaial l|l984l) and iHauschildtl 

we discrctizc the wavelength derivative in the 
SSRTE with a fully implicit method in order to ensure 
stability: 



dXI 
~dX 



Xd. 



x, 



so that the wavelength-discretized SSRTE becomes 



dl\, A//a, — A/_i/a ; _i 



V\i - (xx l +4a A )J Ai . 



ds Xi - A;_i 

If we now define the optical depth scale along the ray 



as 



dr = x + ax 4 



A, 



= Xds, 



A; — Xi-t , 

and introduce the source function, S = Tj/Xi we have 
a\ Xi-i 



dr x 



X Xi - A 



i-i 



'Ai_i 



= I-S. 



With this definition, the formal solution of the SSRTE 
along t he characteristics can b e written in the following 
way (cf. lOlson fc Kunaszlll987L for a derivation of the for- 
mulae) : 



I{n) = J(Ti_i)exp(Ti_i —n) 
S(r)exp(ri_i - 



I(n) = 7 J _ 1 exp(-AT l _ 1 ) + Ali 



(6) 



(7) 



where we have suppressed the index labeling the ray; Tj 
denotes the optical depth along the ray with T\ = and 
rj_i < rj while r is calculated using piecewise linear in- 
terpolation of x along the ray, viz. 



At,:. 



(Xi-i +Xi)\si 



8i\/2. 



(8) 
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Fig. 1. Absorptive continuum (e c = 1). Top part: CMF 
H in arbitrary units (full line: £ — 0, dotted line: ( — 
1), bottom part: relative differences between CMF H (full 
line) and CMF J (dotted line) for the two discretizations. 
The + signs give the location of the actual wavelength 
points used in the computation. 

The source function S (r) along a characteristic is interpo- 
lated by linear or parabolic polynomials so that 



AJj = ctjSj-i + (3 t Si + 7i<Si + i, 



(9) 



Ari = T,-|_i — Ti is the optical depth along the ray from 
point i to point i + 1. The coef ficients a^, /3j, and ji are 
given in lOlson fc Kunaszl l|l987t) . 

For the tangent rays (for example ray 1 in Fig. 1 
of lHa,iischildl h992V the formal solution starts at point 
2 with Ji given as the outer boundary condition and 
proceeds along the ray. The formal solution for a core- 
intersecting ray is split into two parts: (i) integration from 
point 1 to point N, where I\ is given as the outer boundary 
condition and (ii) integration from point N + 2 to point 
27V, where /jv+ i is given as the inner boundary condition. 

2.1.2. Second dl/dX discretization 

The discretization of the 81/ dX derivative in the previous 
section can be deferred. We first rewrite Eq. [3] as 



dh _ ai dXI _ j _ W_ 
dr xi dX xi 



(10) 



with 

Xl = Xi 
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and the definition of the CMF optical depth 
dr = -xids 

along the characteristic. To obtain an expression for the 
formal solution, we rewrite Eq. ^| as 

^ = i-s-s 

dr 
with 

s=*s= r ± 



(11) 

S = —'— (12) 

(13) 

With this we obtain the follo wing expression fo r the 
formal solution (see also Eq. 20 in lHauschildtll 1 992^ 



XI XI 
ai dXI 



h,i = Ii-i,iexp{-Ari-i] 
with the definitions 



SliA + SL 



(14) 



Slid = atijSi-ij + /9i,j5i,j +7i,;-5j+i,i 
and 

5h.i = a^iSi-iA + PijSij 

The index i labels the (spatial) points along a character- 
istic, the index I denotes the waveleng th point. The coef - 
ficie nts an, 0n, and 7» i a re given in IHauschildtl l)l992|) 
and lOlson fe Kunasd l)l987|) . here they are calculated for a 
fixed wavelengths for all points along a characteristic. S is 
a vector of known quantities (the old mean intensities and 
thermal sources). The S contain the effects of the velocity 
field on the formal solution and are given by 

- a,_i,i dXI 
Ji-l.l = —-z 



Xi-l,l 



Sii — — 



a t , dXI 



dX 



Note that the integration of S is performed using linear 
elements in order to allow for a recursive formal solution. 
Parabolic elements can be used but require the solution 
of matrix equations to obtain the formal solution. 

If the velocity field is monotonically increasing or de- 
creasing we use a stable upwind discretization of the wave- 
length derivative. In both cases, the problem becomes an 
initial value problem and can be solved for each wave- 
length once the results of the previous (smaller or longer) 
wavelength points are known. For a monotonically increas- 
ing velocity field this gives 

- a>i-i,i Xi 

-1,1 — ~ J J-i-1,1 



Xi-l,l 

Xi- 



Xi-X 



i-i 



Si. 



Xi 



Xi- 



--fi-lJ-l 



1 

A, 



Ai-i 



-hi 



A; — A 



l-l 



(15) 
(16) 
(17) 
(18) 



Here, and in the following, we use a sorted wavelength grid 
with A/_i < A; < A/+i for monotonically increasing veloc- 
ities (the order of the grid would be reversed for monotoni- 
cally decreasing velocity fields). For non-monotonic veloc- 
ity fields this is no longer the case and the formal solution 
needs to explicitly account for the wavelength couplings 
I Bar on fc HauschildtjEiol . 



With these formulae we can write 5Ii y i in the form 

SliA = &iA (pi-Xdk-1,1 - Pi-lJ-lIi-lA-l) 
(pi,lli,l - Pi,l-\h,l-l) 



where 



Pi- 1,1 = 

Pi-i,l-i = 

Pi,l = 

PU-l = 



X, 



A/ — A;_i Xi-l,l 

Xi-\ ai-i,i 

A/ — A;_i Xi-1,1 

A; aid 

A/ — A;_i Xid 

A;-i aid 



A; - A;_i x 



id 



The formal solution then assumes the form 
(l - j3idPi,i) Ii,i = (&id +exp(- Arj_i ))h-\,i 

— OlidPi-ld-lIi-ld-l 
-0i,lPi,l-lh,l-l + Sh,l 



This equation can be solved recursively along a character- 
istic to calculate the I id for all i and fixed I. The form of 
the wavelength discretization is now second order accurate 
in AA. 

The new form of the formal solution requires only small 
ch anges to the constr uctio n of the A* operator discussed 
in IHauschildtl l)l992|) and lHauschildt. Storzer. fc Baronl 
(1991 - 



We describe the construction of A* for arbitrary band- 
width using the example of a tangential characteristic. The 
intersection points (including the point of tangency) are 
labeled from left to right, the direction in which the formal 
solution proceeds. For convenience, we label the character- 
istic tangent to shell k + 1 as k. Therefore, the character- 
istic k has 2k + 1 points of intersection with discrete shells 
1 . . . k + 1. To compute row j of the discrete A-operator 
(or A-matrix), Ay , we sequentially label the intersection 
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points of the characteristic k with the shell i and defi 
auxiliary quantities A& and as follows: 



mc 



\ k 



\ k 

A i-i 



"*3 **"" y 
for i < j — 1 



= 



for i = j — 1 



-l 



for i = j 




**p*_ ltI + expt-AT^U) +/i 



- (1-/3^ 



for i = j + 1 

for j + 1 < i < k + 1 
1 



*3+l 



I Pi. I 

,fc„fe 



Pi-x.1 +cxp(-Arfl 1 ) 
For the calculation of X k j, we obtain: 



\ k 



\ k 



\ k 



for i = k + 2 

[a J fc Pti : z+ ex P(- AT f-i) 
for fc + 2 < i < fc + j + 2 



+ ex p(- Ar ti) 

for j = fc + j + 2 



■ /f Pi,* 



-1 



\-l,3 



a k -" k 



Pti,i+ ex P(- Ar f-i! 



for i = fc + j + 3 

x -1 



I \ k 

1 \-lJ 



&Pti,i +exp(-Arf_ x ; 



for i = fc + j + 4 



for fc 



% - (l - 



j + 5 < i < 2k + 1 
1 



I \ k 

1 A i-i,j 



-exp(-Arf_ i ; 



Using the A^ and A^ , we can now write the A- Matrix as 



A., 



{0 



where w k ^ are the angular quadrature weights, {1} is the 
set {i < k + 1} and {Z'} is the set {i > k + 1}. This expres- 
sion gives the full A-matrix, it can easily be specialized to 
compute only certain bands of the A-matrix. In that case, 



not all of the X k j and X k j have to be computed, reducing 
the CPU time from that required for the computation of 
the full A-matrix. 



2.1.3. Stability 

In the second discretization scheme, the wavelength 
derivative contains an explicit term, which is required to 
derive a recursive method with second order accuracy. We 
can stabilize the discretion via a method similar to the 
stand ard Crank-Nicholson scheme ijAbramowitz k, Stegunl 
119721) . To combine the two discretization schemes de- 
scribed above we i ntroduce a factor C € [0, 1] so that for 
the first scheme in lHauschildtl l)l992|) we have C = 1 an d 
for the second discretization we have C = 0. With this the 
optical depth scale along the ray dr and x become 



dr = x + a x 4 + C 



' A, - A 



1-1 



and S is given by 



s + c- 



X A; - A 



1-1 



whereas S becomes 
ai dXI 



S 



-(i-C) 



X dX 



With this we can arbitrarily combine the two dl/dX dis- 
cretizations by varying £ over the interval [0, 1]. This al- 
lows us to combine them to optimize accuracy of the over- 
all solution of the SSRTE. 

2.1.4. Discussion 

It should be expected that the two discretization schemes 
introduced above will show different numerical behavior. 
The C — 1 scheme handles the Ix [ and /ai_i parts dif- 
ferently since the (known) is treated like a source 
term whereas the Ix t term leads to a changed definition 
of dr in the CMF along the ray. Therefore, the I\ t is ba- 
sically assumed to vary exponentially across a character- 
istic. This will introduce more diffusive behavior for large 
optical depths and strong scattering (where the intrinsic 
assumptions for Ix, and will differ most). 

This is not the case for the £ = discretization which 
assumes that the Xdl/dX itself can be fitted to a first order 
polynomial and treated as a source term. The discretiza- 
tion step is delayed until the integration of the source 
term. Although the discretization is implicit by design, 
the behavior of 81/ dX is a priori unknown and the sec- 
ond order nature of the discretization leads to an explicit 
term. Under extreme conditions, e.g., a strong line in an 
optically thin rapidly expanding gas, it is possible that the 
numerical integration is dominated by the "known" term 
due to large changes in 81/ dX. This leads to numerical in- 
stability. We therefore expect that the two discretizations 
will show different behavior in numerical tests. 
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3. Tests 

In order to compare and test the two different discretiza- 
tion schemes, we have devised a simple test model. We 
consider the case of a 2-level atom with background con- 
tinuum on a radial grid (assuming spherical symmetry). 
The background continuum is assumed to be grey with 
adjustable thermalization fraction e c = n c /Xc with Xc = 
K c + a c . The line of the 2-level atom is parameterized by 
its strength relative to the continuum Xl/Xc and its ther- 
malization fraction e; — ki/xi- F° r the intrinsic line profile 
we use a Gaussian profile with a AAd of 30kms _1 . For 
the radial structure, we set %c oc r~ 2 and use a radial grid 
with i?out = 101 and i?; n = 1 (so that the extension of the 
test atmosphere is 100) and a logarithmic optical depth 
grid from r = 10 4 to 10 -6 in the continuum to fix the 
scale of the opacities. The velocity law is linear (homol- 
ogous expansion) with a prescribed maximum velocity of 
1000 km s -1 . The computations are performed on a wave- 
length grid that uses a stepsize of 0.1 line width inside 
the line and 5 times the line width outside the line if not 
specified in detail for some of the test calculations. 



3.1. Continuum tests 

As a first test we verified that there are no differences 
between the solutions for £ = and £ = 1 for zero ex- 
pansion velocity. The next test is for a pure continuum 
case (i.e., line strength of zero). This is actually a difficult 
test for the method as the flat continuum should be repro- 
duced by the algorithm. The results for purely absorptive 
(e c = 1) and scattering dominated (e c = 10~ 2 ) continua 
are shown in Figs^ and [21 respectively. In this case the 
CMF H (the first Eddington moment of the intensity in 
the CMF) should be constant. From Fig. ^ we see that 
both discretization schemes deliver nearly identical results 
with differences below 0.3% in the case of an absorption 
dominated continuum. For the case of a scattering dom- 
inated continuum the differences are significantly larger, 
nearly 6%, cf. Figure [5J From the top part of the figure it 
is obvious that the differences start to increase just around 
the region of resolution change in wavelength at the (zero 
strength) line. That scattering increases the diffusive ef- 
fect of the C — 1 wavelength discretization as shown by 
comparing Figs. ^ and [5] This effect becomes larger and 
more pronounced if the wavelength resolution is lower, 
which is shown in Figs. [3] and 0] In this test, we have 
decreased the wavelength resolution by roughly a factor 
of 10. The case of the absorption dominated continuum 
shows similar behavior as the test with higher wavelength 
resolution with the relative differences increasing slightly. 
However; the scattering dominated case produces much 
larger relative differences of up to 13% between the two 
discretizations. Whereas the C — discretization produces 
the expected flat continuum (once the initial conditions in 
wavelength are "forgotten" ) , the £ = 1 discretization pro- 
duces a pseudo line feature just due to the change in res- 
olution. This illustrates the higher accuracy of the second 




Fig. 2. Scattering dominated continuum (e c = 10~ 2 ). Top 
part: CMF H in arbitrary units (full line: ( = 0, dot- 
ted line: £ = 1), bottom part: relative differences between 
CMF H (full line) and CMF J (dotted line) for the two 
discretizations. The + signs give the location of the actual 
wavelength points used in the computation. 



order scheme, which since it is properly centered produces 
more accurate results on an irregular grid. 

The case of a constant step size wavelength grid is 
shown in Fig. [S] The step size in wavelength is 300 km s -1 . 
The £ = discretization produces a perfectly flat contin- 
uum after the effect of the (grey) initial condition dies 
away. In contrast, the £ = 1 discretization produces con- 
siderable numerical diffusion in this test, on the average 
the solution deviates by more than 10% from the £ = 
solution. In addition, the continuum is not flat but shows 
a noticeable increase in CMF H (and J) toward the red. 

It is clear that the £ = discretization gives signif- 
icantly better results than the £ = 1 approach for the 
continuum tests. 

3.2. Line tests 

In contrast to the previous set of tests, we now introduce a 
strong line. We set the line strength Xl/Xc — 100 with e; = 
10 -4 to simulate a strong, scattering dominated line of a 
2-level atom against a background continuum. Figs. Eland 
Ejshow the results for absorption and scattering dominated 
background continua. In both cases the relative differences 
in the continua are comparable to the previous tests. The 
differences in the lines are actually comparatively small in 
both cases, significantly smaller than might be expected 
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Fig. 3. Absorptive continuum (e c = 1). Top part: CMF 
H in arbitrary units (full line: £ = 0, dotted line: ( = 
1), bottom part: relative differences between CMF H (full 
line) and CMF J (dotted line) for the two discretizations. 
The + signs give the location of the actual wavelength 
points used in the computation. 



Fig. 5. Scattering dominated continuum (e c = 1CP 2 ). Top 
part: CMF H in arbitrary units (full line: £ = 0, dot- 
ted line: £ = 1), bottom part: relative differences between 
CMF H (full line) and CMF J (dotted line) for the two 
discretizations. In this calculation equidistant wavelength 
grid with a spacing of 300 km s^ 1 was used. 



-'0000 -5D0D 



scattering con 




Fig. 4. Scattering dominated continuum (e c = 10~ 2 ). Top 
part: CMF H in arbitrary units (full line: C = 0, dot- 
ted line: £ — 1), bottom part: relative differences between 
CMF H (full line) and CMF J (dotted line) for the two 
discretizations. The + signs give the location of the actual 
wavelength points used in the computation. 



from the continuum test results considering the fact that 
the line is strongly scattering dominated. This is most 
likely caused by the smoothing effect of the rapidly varying 
opacity across the line. The effect of continuum scattering 
is to considerably widen the line due to line photons being 
scattered by the continuum. The more diffusive £ = 1 
discretization scheme produces an emission feature that is 
about 2-3% stronger than the C = discretization relative 
to the continuum. It is interesting to note the strong effect 
of even a weak (compared to the line itself) background 
continuum on the shape of the line and the differences 
between the two discretizations. 

For the final test we show in Fig. he results of a 
calculation for £ = 0.5. Compared to £ = I this calculation 
shows much smaller differences from the £ = case in the 
continuum. From the top portion of the figure it is clear 
that in the line itself the £ = 0.5 model falls roughly half 
way in between the £ — 1 and £ — cases. 

In all cases the £ = discretization gives numerically 
better results. However, there are situations where the £ = 
discretization is numerically unstable, which we discuss 
in the following section. 

3.3. Numerical instability for £ = 

In certain conditions, it is possible for the £ = discretiza- 
tion to become numerically unstable. This is illustrated 
in the test case shown in Fig. |H| The two test lines are 
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absorptive continaam 




-1.0 I i I i i i I i i i I i i i I i i i I , 

-4000 -2000 2000 4000 

distance from line centen [km/.] 

Fig. 6. Strong (xi/Xc = 100) scattering dominated (e; = 
10~ 4 ) line with an absorptive background continuum (e c = 
1). Top part: CMF H in arbitrary units (full line: C 
0. dotted line: £ = 1), bottom part: relative differences 
between CMF H (full line) and CMF J (dotted line) for 
the two discretizations. 

in complete LTE (S = B), the background continuum is 
scattering dominated. The structure is taken from typi- 
cal SN la structure, but we have completely ignored the 
strong line blanketing inherent in SNe la by including 
line opacity from only Ca II. This procedure is actually 
quite useful for line identification in detailed model calcu- 
lations. The velocity field is homologous (v oc r) with a 
maximum speed of 30, 000 km s -1 , the inner boundary is 
at v = 2,000kms _1 . We use this model rather than the 
simpler more academic test cases for two reasons: 1) this 
is the model where we actually discovered the instability 
for C = 0, and 2) we have been unable to find a simple 
test case that so strongly reproduces the instabilities. The 
continuum is optically thin in absorption (r a b s ~ 4 x 10 -4 ) 
and reaches a scattering optical depth T scat t « 1.5. We use 
a Planck-function as the inner boundary condition for the 
intensities for simplicity, physically it is better to employ 
a nebular [symmetry] boundary condition which is shown 
in Figure We see that the instability is evident (for the 
C = case) even though the line has gone into emission 
(note that the different scales between Fig. [S] and Fig. [5] 
give the impression that the instability is somewhat sup- 
pressed in the nebular case, but it is simply due to the 
larger range required to plot the strong emission features). 
The unstable wiggles are not associated with the grid and 
varying the wavelength resolution does not significantly af- 
fect the output spectrum. While we believe from our tests 



scattering cantinuam 
4xi o 9 r 1 i 1 1 1 i 1 1 ^ 1 1 n~ 




- 



-10 I , I , , , I , , , I , , , I , , , I , 

-4000 -2000 2000 4000 

distance Irom line cente, (km/s) 

Fig. 7. Strong (xi/Xc — 100) scattering dominated (e; = 
10 ) line with a scattering dominated continuum (e c = 
10~ 2 ). Top part: CMF H in arbitrary units (full line: C = 0, 
dotted line: £ = 1, dashed line: £ = 0.5), bottom part: rela- 
tive differences between CMF H for the two discretizations 
with £ = and £ = 1. In addition, the dotted line shows 
the relative differences between calculations with £ = 0.5 
and C = 0. 

that the instability is produced predominantly in optically 
thin atmospheres, we have not been able to ascertain the 
precise conditions that trigger the instability. 

The instability in the £ = discretization is obvious 
starting at the rest wavelength of the line. Even in the 
line trough the oscillations created by the instability are 
apparent. The overall line shapes are distorted, even in the 
red emission feature. The instability vanishes for £ ^ as 
shown in the Fig. |S| The differences between the spectra 
for £ = 0.5 and 1.0 are very small, the £ = 0.1 case shows 
a few percent difference compared to £ = 1. It is clear that 
the instability of the £ = case can easily be avoided by 
using C > due to the stabilizing features of the £ = 1 
discretization. In our test case even £ = 0.1 is sufficient to 
prevent the instability. 

4. Conclusions 

We have studied the discretization of the -Jr term in 
the co-moving frame radiative transfer equation. We have 
found that we can write a second order discretization 
scheme, which is in general more accurate than our previ- 
ous method, but it is possible for this scheme to become 
unstable. We have constructed a hybrid Crank-Nicholson 
like scheme, which is unconditionally stable. We have been 
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C\T Wavelength [A] 

Fig. 8. Example to illustrate the possibility of numerical 
instabilities in the £ = discretization. The structure is 
taken from an Type la supernova calculation, only Ca II 
lines are included in the opacity and the lines are assumed 
to be in complete LTE. 

io i6 F ' ' ' ' I ' ' ' ' I ' ' ' ' I ' 1 ' ' E 



( = 0.0 
'=[ I 

}= '.o 




Fig. 9. The same structure as was used in Fig. 03 but with 
nebular (symmetry) boundary conditions. 

unable to identify the exact conditions which lead to the 
unstable behavior, but we have shown that even with a 
small admixture of the unconditionally stable scheme, sta- 
bility is recovered. We recommend using a small value of 
£ = 0.1, which in all our tests leads to recovered stability 
and is more accurate than a higher value of £. However, 
until we determine the exact conditions that trigger the 
instability the cautious user should vary £ and determine 
the sensitivity of the results to its value. 
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